function [fitdat]= lambdas(tsdat,stepdat,coren,nsteps,titlestr)
%function [fitdat]= lambdas(tsdat,stepdat,coren,nsteps,titlestr)
%
% modified from lambda.m - originally by Steve Suttles
% modified by Patrick Dickhudt Dec. 2010
% - removed evals
% - simplify plotting
% - add expfit subfunction for fitting, eliminating need for curve fitting
% toolbox
% - change x values in fit to step elapsed time rather than expt. elapsed
% time. improves robustness of fitting and makes E0 from fit meaningful
%
% Changes
% 13-Apr-2011 - change calculation of ylims for plt_lambda subfunction to
% exclude nans - PJD

for i1=1:nsteps
    
    [efit,r2,tmin,tmax] = plt_lambda(tsdat,stepdat,coren,i1,titlestr);
    
    close
    
    fitdat(i1).stepn = i1;
    fitdat(i1).E0 = efit(1);
    fitdat(i1).lambda = efit(2);
    fitdat(i1).r2 = r2;
    fitdat(i1).tmin = tmin;
    fitdat(i1).tmax = tmax;
    
    fitdat(i1).fields = {'fitdat = results of exponential fit for each step';
         'stepn = step number';
         'E0 = initial erosion rate from fit (kg m-2 s-1), ~ E(tmin)';
         'lambda = exponential decay constant for each step (s-1)';
         'r2 = r-squared statistic'
         'tmin = beginning elapsed time used in fit';
         'tmax = ending elapsed time used in fit'};

end


                 
%save([coren '_lambda'],'fitdat')

%--------------------------------------------------------------------------

function [efit,r2,tmin,tmax]=plt_lambda(tsdat,stepdat,core,step,titlestr)
%function [cf_,gof_,tmin,tmax]=plt_lambda(tsdat,stepdat,num,exp,core,step,titlestr)
%m-file to plot ea step and find lambda values doing exponential fits

% get erosion rate data for current step
y_E = tsdat.E(stepdat.idx.(['step' num2str(step)])); 

% get time data for current step
x_etime = tsdat.etime(stepdat.idx.(['step' num2str(step)])) - min(tsdat.etime(stepdat.idx.(['step' num2str(step)])));

tmin=x_etime(1);
tmax=x_etime(length(x_etime));

figure

cont = 'n';

while cont == 'n'
    [efit,model,r2]=fit_lambda(x_etime,y_E,tmin,tmax);
    [sse,FittedCurve] = model(efit);

    figure(gcf)
    clf;
    
    set(gcf,'name',['Lambda fit  Core ' titlestr '  Step' num2str(step)])
    iplt=find(x_etime >= tmin & x_etime <= tmax);
    size(iplt);
    p = plot(x_etime,y_E,'-b',x_etime(iplt),FittedCurve,'-r');
    set(p,'linewidth',1.2)
    set(gca,'yscale','log')
    ylow = 10^(floor(log10(min(y_E(y_E > 0))))-1);
    yhigh = 10^(ceil(log10(max(y_E(y_E > 0))))+1);
    set(gca,'ylim',[ylow yhigh])
    title(['Lambda fit  Core ' titlestr '  Step' num2str(step)])
    xlabel('Step elapsed time, [s]')
    ylabel('E, [kg ^. m^{-2} ^. s^{-1}]')
    eqnstr =['E = ' num2str(efit(1),'%5.3e') ' ^. e ^{-' num2str(efit(2),'%8.5f') ' * t}'];
    xl=[x_etime(1) x_etime(length(x_etime))];
    ylim=get(gca,'ylim');
    xlim(xl)
    ytp1 = 10^(log10(ylim(2)) - (log10(ylim(2)) - log10(ylim(1)))*.1);
    ytp2 = 10^(log10(ylim(2)) - (log10(ylim(2)) - log10(ylim(1)))*.2);
    text(xl(1)+0.1*(xl(2)-xl(1)),ytp1,eqnstr)
    r2str = ['r^2 = ' num2str(r2,'%5.3f')];
    text(xl(1)+0.1*(xl(2)-xl(1)),ytp2,r2str)
    lbstr = ['t_{min} = ' num2str(round(tmin))];
    ubstr = ['t_{max} = ' num2str(round(tmax))];
    text(xl(1)+0.75*(xl(2)-xl(1)),ytp1,lbstr)
    text(xl(1)+0.75*(xl(2)-xl(1)),ytp2,ubstr)

    disp(['Step end time = ' num2str(xl(2))])
    disp(' ')
    
    figure(gcf)
    cont = input('accept this fit, [y/n]: ','s')

    if cont == 'n'
        disp(' ')
        tmin=input('enter new lower limit of fit :')
        disp(' ')
        tmax=input('enter new upper limit of fit :')
    end
end

saveas(gcf,[core '_step' num2str(step) '_lambdafit'],'fig')
            
end

%--------------------------------------------------------------------------

function [efit,model,r2]=fit_lambda(x_etime,y_E,xmin,xmax)
%FIT_LAMBDA    Create plot of datasets and fits
%   FIT_LAMBDA(X_ETIME,Y_E)
%   Creates a plot, similar to the plot in the main curve fitting
%   window, using the data that you provide as input.  You can
%   apply this function to the same data you used with cftool
%   or with different data.  You may want to edit the function to
%   customize the code and this help message.
%
%   Number of datasets:  1
%   Number of fits:  1

 
% Data from dataset "y_E vs. x_etime":
%    X = x_etime:
%    Y = y_E:


x_etime = x_etime(:);
y_E = y_E(:);

% indices of data to include in fit
ftinds = xmin<=x_etime & x_etime<=xmax & ~(isnan(x_etime) | isnan(y_E));

% estimate starting points for fit using data
stp1 = y_E(min(find(ftinds)));
stp2 = -((mean(log(y_E(ftinds))) - log(stp1)))./ floor(mean(find(ftinds)));

% fit exponential function to time vs TSS data
% Fit this model using new data
if sum(ftinds)<2  %% too many points excluded
   error('Not enough data left to fit ''%s'' after applying exclusion rule ''%s''.','fit 1','excl')
else
    [efit,model] = expfit(x_etime(ftinds),y_E(ftinds),[stp1 stp2]);
end

[sse, FittedCurve] = model(efit);
sst = sum((y_E(ftinds) - mean(y_E(ftinds))).^2);
r2 = 1-sse./sst;

function [estimates, model] = expfit(xdata, ydata,start_point)
    % exponential fit function using fminsearch
    % copied from Matlab help: Curve Fitting via Optimization
    % renamed, in help called fitcurvedemo, conveniently fittype was exactly
    % right for this purpose
    %
    % added by Patrick Dickhudt 2-Dec-2010

    % Call fminsearch with a random starting point.
    %start_point = [.005 .0005];
    model = @expfun;
    options = optimset('TolFun',1E-6,'TolX',1E-6,'MaxFunEvals',1000,'MaxIter',1000);
    [estimates,fval,exitflag] = fminsearch(model,start_point,options);
    
    % expfun accepts curve parameters as inputs, and outputs sse,
    % the sum of squares error for A*exp(-lambda*xdata)-ydata,
    % and the FittedCurve. FMINSEARCH only needs sse, but we want
    % to plot the FittedCurve at the end.
    
function [sse, FittedCurve] = expfun(params)
        A = params(1);
        lambda = params(2);
        FittedCurve = A .* exp(-lambda * xdata);
        ErrorVector = FittedCurve - ydata;
        sse = sum(ErrorVector .^ 2);
end
end


end
end
%--------------------------------------------------------------------------